Exponentiated Weibull distribution (exponweib)

The exponentiated Weibull distribution (SciPy: scipy.stats.exponweib) is a flexible continuous distribution on \([0,\infty)\) (in its standard form) obtained by raising the Weibull CDF to a positive power.

It is popular in reliability and survival analysis because—with two shape parameters—it can represent a wide range of hazard-rate behaviors (increasing, decreasing, unimodal, and bathtub-like).

---

Learning goals

• Write down the PDF/CDF/PPF and connect them to the Weibull distribution.

• Understand how the parameters change shape and hazard behavior.

• Derive moments using a clean change-of-variables + series expansion.

• Implement inverse-CDF sampling using only NumPy.

• Use SciPy (exponweib) for evaluation, simulation, and MLE fitting.

import numpy as np

import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

import scipy
from scipy import integrate, optimize, special
from scipy.stats import exponweib as exponweib_dist
from scipy.stats import weibull_min, chi2, kstest

# Plotly rendering (CKC convention)
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

# Reproducibility
rng = np.random.default_rng(42)
np.set_printoptions(precision=6, suppress=True)

# Record versions for reproducibility (useful when numerical details matter)
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}
VERSIONS

{'numpy': '1.26.2', 'scipy': '1.15.0', 'plotly': '6.5.2'}

1) Title & Classification

• Name: exponweib (Exponentiated Weibull; SciPy: scipy.stats.exponweib)

• Type: Continuous

• Support (standard form): \(x \in [0,\infty)\)

• Parameter space (standard form): shape parameters \(a>0\), \(c>0\)

• SciPy location/scale: loc \in \mathbb{R}, scale > 0 with

  \[X = \text{loc} + \text{scale}\,Y, \qquad Y \sim \mathrm{ExponWeib}(a,c).\]

Unless stated otherwise, we work with the standard form (loc=0, scale=1).

Notation note: many texts write the parameters as \((\alpha, k, \lambda)\), where \(\alpha\) is the exponentiation parameter, \(k\) is the Weibull shape, and \(\lambda\) is the scale. SciPy uses \((a, c)\) for the two shape parameters and scale for \(\lambda\).

2) Intuition & Motivation

What it models

Start with a Weibull random variable \(W\) (shape \(c\), scale \(\lambda\)):

\[F_W(x) = 1 - \exp\left(-(x/\lambda)^c\right), \qquad x \ge 0.\]

The exponentiated Weibull distribution raises this CDF to a power \(a>0\):

\[F_X(x) = \bigl(F_W(x)\bigr)^a = \left[1 - \exp\left(-(x/\lambda)^c\right)\right]^a.\]

So it is literally the Weibull CDF, exponentiated.

A useful interpretation when \(a\) is an integer (\(a=m\in\mathbb{N}\)):

• If \(W_1,\dots,W_m\) are i.i.d. Weibull\((c,\lambda)\), then $\(\max\{W_1,\dots,W_m\} \sim \mathrm{ExponWeib}(a=m, c, \lambda).\)$

This “maximum-of-\(m\) Weibulls” story helps build intuition: larger \(a\) pushes mass to the right (stochastically larger lifetimes).

Typical real-world use cases

• Reliability engineering: modeling lifetimes with non-monotone hazard (e.g. bathtub-shaped failure rates).

• Survival analysis: flexible parametric baseline hazard.

• Hydrology / environmental extremes: positive-valued quantities with skew and flexible tails.

• Manufacturing / materials: time-to-failure with early-life and wear-out effects.

Relations to other distributions

• \(a=1\) gives the standard Weibull distribution (scipy.stats.weibull_min).

• \(c=1\) gives the exponentiated exponential distribution.

• \(a=1\) and \(c=1\) gives the exponential distribution.

• Location/scale transform recovers different units and origins: \(X=\text{loc}+\text{scale}\,Y\).

3) Formal Definition

We present the standardized form first (loc=0, scale=1).

CDF

For \(x \ge 0\), \(a>0\), \(c>0\):

\[ F(x; a,c) = \left[1-\exp\left(-x^c\right)\right]^a. \]

PDF

For \(x>0\):

\[\begin{split} \begin{aligned} f(x; a,c) &= \frac{d}{dx}F(x;a,c) \\ &= a\,c\,x^{c-1}\,\exp(-x^c)\,\left[1-\exp(-x^c)\right]^{a-1}. \end{aligned} \end{split}\]

Quantile function (PPF)

For \(q\in(0,1)\):

\[ F^{-1}(q;a,c) = \left[-\log\left(1-q^{1/a}\right)\right]^{1/c}. \]

This inverse-CDF is the key to fast sampling.

Location / scale

With loc and scale>0, define \(z = (x-\text{loc})/\text{scale}\). Then for \(x\ge \text{loc}\):

\[ F_X(x) = F(z;a,c), \qquad f_X(x) = \frac{1}{\text{scale}}\,f(z;a,c). \]

Hazard function (survival analysis)

The hazard rate is

\[h(x) = \frac{f(x)}{1-F(x)}.\]

For the exponentiated Weibull, the hazard can be increasing, decreasing, unimodal, or bathtub-shaped depending on \((a,c)\).

def exponweib_cdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """CDF of the exponentiated Weibull (NumPy implementation).

    Matches SciPy's `scipy.stats.exponweib.cdf` for the same (a, c, loc, scale).
    """
    if a <= 0 or c <= 0 or scale <= 0:
        raise ValueError("Require a>0, c>0, scale>0")

    x = np.asarray(x, dtype=float)
    z = (x - loc) / scale

    out = np.zeros_like(z, dtype=float)
    mask = z >= 0
    zm = z[mask]

    # exm1c = 1 - exp(-zm**c), computed stably for small zm**c
    exm1c = -np.expm1(-(zm**c))
    out[mask] = exm1c**a
    return out


def exponweib_logpdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Log-PDF of the exponentiated Weibull (stable for small/large x)."""
    if a <= 0 or c <= 0 or scale <= 0:
        raise ValueError("Require a>0, c>0, scale>0")

    x = np.asarray(x, dtype=float)
    z = (x - loc) / scale

    out = np.full_like(z, -np.inf, dtype=float)
    mask = z > 0
    zm = z[mask]

    neg_zc = -(zm**c)
    exm1c = -np.expm1(neg_zc)  # = 1 - exp(-zm**c) in (0,1)

    logp = (
        np.log(a)
        + np.log(c)
        - np.log(scale)
        + (a - 1.0) * np.log(exm1c)
        + neg_zc
        + (c - 1.0) * np.log(zm)
    )
    out[mask] = logp
    return out


def exponweib_pdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """PDF computed from the log-PDF."""
    return np.exp(exponweib_logpdf(x, a, c, loc=loc, scale=scale))


def exponweib_ppf(q: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Percent-point function (inverse CDF)."""
    if a <= 0 or c <= 0 or scale <= 0:
        raise ValueError("Require a>0, c>0, scale>0")

    q = np.asarray(q, dtype=float)
    if np.any((q < 0) | (q > 1)):
        raise ValueError("q must be in [0,1]")

    out = np.empty_like(q, dtype=float)
    out[q == 0] = loc
    out[q == 1] = np.inf

    mask = (q > 0) & (q < 1)
    qm = q[mask]
    out[mask] = loc + scale * (-np.log1p(-(qm ** (1.0 / a)))) ** (1.0 / c)
    return out


def exponweib_rvs_numpy(
    a: float,
    c: float,
    *,
    loc: float = 0.0,
    scale: float = 1.0,
    size=1,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Draw samples from ExponWeib(a,c) using **only NumPy**.

    Inverse transform sampling:
      U ~ Uniform(0,1)
      X = loc + scale * (-log(1 - U^{1/a}))^{1/c}
    """
    if a <= 0 or c <= 0 or scale <= 0:
        raise ValueError("Require a>0, c>0, scale>0")

    rng = np.random.default_rng() if rng is None else rng
    u = rng.uniform(size=size)

    # Avoid log(0) at the endpoints.
    eps = np.finfo(float).eps
    u = np.clip(u, eps, 1.0 - eps)

    return loc + scale * (-np.log1p(-(u ** (1.0 / a)))) ** (1.0 / c)


def exponweib_sf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Survival function S(x)=1-CDF(x), computed stably via log1p."""
    cdf = exponweib_cdf(x, a, c, loc=loc, scale=scale)
    return np.exp(np.log1p(-cdf))


def exponweib_hazard(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Hazard h(x)=f(x)/S(x), computed in log-space for stability."""
    x = np.asarray(x, dtype=float)
    logf = exponweib_logpdf(x, a, c, loc=loc, scale=scale)
    sf = exponweib_sf(x, a, c, loc=loc, scale=scale)
    out = np.zeros_like(x, dtype=float)

    mask = sf > 0
    out[mask] = np.exp(logf[mask] - np.log(sf[mask]))
    out[~mask] = np.inf
    return out

# Sanity check: our formulas match SciPy (prefer logpdf to avoid overflow near 0).

a, c, loc, scale = 1.7, 1.3, 0.4, 2.2
x = np.logspace(-5, 2, 40) * scale + loc
qs = np.linspace(0.01, 0.99, 9)

assert np.allclose(exponweib_logpdf(x, a, c, loc=loc, scale=scale), exponweib_dist.logpdf(x, a, c, loc=loc, scale=scale))
assert np.allclose(exponweib_cdf(x, a, c, loc=loc, scale=scale), exponweib_dist.cdf(x, a, c, loc=loc, scale=scale))
assert np.allclose(exponweib_ppf(qs, a, c, loc=loc, scale=scale), exponweib_dist.ppf(qs, a, c, loc=loc, scale=scale))

# Check that inverse-CDF sampling produces the right distribution (quick KS test).
x_samp = exponweib_rvs_numpy(a, c, loc=loc, scale=scale, size=3_000, rng=rng)
D, p = kstest(x_samp, lambda t: exponweib_dist.cdf(t, a, c, loc=loc, scale=scale))
print(f"KS statistic={D:.4f}, p-value={p:.4f}")

KS statistic=0.0092, p-value=0.9576

4) Moments & Properties

Raw moments

Assume the standard form with loc=0 and scale=\lambda.

Using the substitution \(y = (x/\lambda)^c\) (so \(x=\lambda y^{1/c}\)), the exponentiated Weibull density simplifies nicely and we obtain a series expression for the raw moments:

\[ \mathbb{E}[X^r] = \lambda^r\,a\,\Gamma\left(1+\frac{r}{c}\right) \sum_{j=0}^{\infty} (-1)^j\,\binom{a-1}{j}\,\frac{1}{(j+1)^{1+r/c}}. \]

• For integer \(a=m\in\mathbb{N}\), the binomial series terminates at \(j=m-1\).

• For non-integer \(a\), the series is infinite but often converges quickly.

From the raw moments we get:

• mean: \(\mu = \mathbb{E}[X]\)

• variance: \(\sigma^2 = \mathbb{E}[X^2] - \mu^2\)

• skewness and kurtosis from the central moments \(\mu_3,\mu_4\).

MGF / characteristic function

There is no simple elementary closed form for the MGF in general. Existence depends on the tail:

• If \(c>1\), the tail is super-exponential and \(M_X(t)=\mathbb{E}[e^{tX}]\) exists for all real \(t\).

• If \(c=1\) (exponentiated exponential tail), \(M_X(t)\) exists for \(t < 1/\lambda\).

• If \(0<c<1\), \(M_X(t)\) diverges for any \(t>0\) (but the Laplace transform exists for \(t\le 0\)).

The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) always exists and can be computed numerically.

Entropy

The differential entropy is

\[H(X) = -\mathbb{E}[\log f(X)].\]

SciPy provides exponweib.entropy(...) (numerical integration). A Monte Carlo estimate using -mean(logpdf(samples)) is also useful for sanity checks.

def exponweib_raw_moment_series(
    r: float,
    a: float,
    c: float,
    *,
    scale: float = 1.0,
    max_terms: int = 50_000,
    tol: float = 1e-12,
) -> float:
    """Compute E[X^r] via the binomial/Gamma series (loc=0).

    Uses the identity:
      E[X^r] = scale^r * a * Gamma(1 + r/c) * sum_{j>=0} (-1)^j * binom(a-1, j) / (j+1)^{1 + r/c}

    Notes
    -----
    - Requires a>0, c>0, scale>0.
    - For integer a, the sum is finite (j=0,...,a-1).
    - For non-integer a, the sum is infinite; we truncate by `tol`.
    """
    if a <= 0 or c <= 0 or scale <= 0:
        raise ValueError("Require a>0, c>0, scale>0")
    if r <= -c:
        raise ValueError("Moment diverges for r <= -c in the standard form")

    p = 1.0 + r / c
    prefactor = (scale**r) * a * special.gamma(1.0 + r / c)

    s = 0.0
    for j in range(max_terms):
        coeff = special.binom(a - 1.0, j)
        term = ((-1.0) ** j) * coeff / (j + 1.0) ** p
        s_new = s + term

        # Stop once the marginal contribution becomes negligible.
        if j > 0 and abs(term) < tol * max(1.0, abs(s_new)):
            s = s_new
            break
        s = s_new

    return float(prefactor * s)


def exponweib_mvsk_from_raw_moments(a: float, c: float, *, scale: float = 1.0):
    """Return mean, variance, skewness, excess kurtosis via raw moments."""
    m1 = exponweib_raw_moment_series(1.0, a, c, scale=scale)
    m2 = exponweib_raw_moment_series(2.0, a, c, scale=scale)
    m3 = exponweib_raw_moment_series(3.0, a, c, scale=scale)
    m4 = exponweib_raw_moment_series(4.0, a, c, scale=scale)

    mean = m1
    var = m2 - m1**2

    mu3 = m3 - 3 * m2 * mean + 2 * mean**3
    mu4 = m4 - 4 * m3 * mean + 6 * m2 * mean**2 - 3 * mean**4

    skew = mu3 / (var ** 1.5)
    kurt_excess = mu4 / (var**2) - 3.0

    return mean, var, skew, kurt_excess

# Compare the series-based moments to SciPy's numerical stats.

a, c, scale = 1.6, 1.2, 2.0

mean_s, var_s, skew_s, kurt_s = exponweib_mvsk_from_raw_moments(a, c, scale=scale)
mean_sp, var_sp, skew_sp, kurt_sp = exponweib_dist.stats(a, c, scale=scale, moments="mvsk")

print("Series moments:")
print("  mean=", mean_s)
print("  var =", var_s)
print("  skew=", skew_s)
print("  kurt(excess)=", kurt_s)

print("\nSciPy stats (numerical):")
print("  mean=", float(mean_sp))
print("  var =", float(var_sp))
print("  skew=", float(skew_sp))
print("  kurt(excess)=", float(kurt_sp))

Series moments:
  mean= 2.4294714656111736
  var = 2.7333895473411074
  skew= 1.3022259281402193
  kurt(excess)= 2.474111034724598

SciPy stats (numerical):
  mean= 2.429471462970749
  var = 2.7333895579577634
  skew= 1.30222591791057
  kurt(excess)= 2.4741110129320223

# Entropy: SciPy (numerical integration) vs Monte Carlo.

a, c, scale = 1.6, 1.2, 2.0
rv = exponweib_dist(a, c, scale=scale)

H_scipy = rv.entropy()

x = rv.rvs(size=60_000, random_state=rng)
H_mc = -np.mean(rv.logpdf(x))

print(f"Entropy (SciPy)      = {H_scipy:.6f}")
print(f"Entropy (Monte Carlo)= {H_mc:.6f}")

Entropy (SciPy)      = 1.756397
Entropy (Monte Carlo)= 1.757041

# Characteristic function (CF) example via Monte Carlo

a, c, scale = 1.2, 0.9, 1.5
rv = exponweib_dist(a, c, scale=scale)

x = rv.rvs(size=120_000, random_state=rng)

def cf_mc(t: float) -> complex:
    return complex(np.mean(np.exp(1j * t * x)))

for t in [0.5, 1.0, 2.0, 4.0]:
    print(f"t={t:>3}: phi(t)≈ {cf_mc(t)}")

t=0.5: phi(t)≈ (0.5638539622591024+0.489223667427389j)
t=1.0: phi(t)≈ (0.2521834895613472+0.4322972648554105j)
t=2.0: phi(t)≈ (0.07024447912956258+0.2711293646917929j)
t=4.0: phi(t)≈ (0.01435532670293679+0.14062825155057773j)

5) Parameter Interpretation

Think of the parameters as controlling (i) the underlying Weibull shape and (ii) how strongly we “exponentiate” its CDF.

Shape parameter \(c\) (Weibull shape)

• \(c<1\): heavy mass near 0, decreasing hazard (infant mortality)

• \(c=1\): constant hazard (exponential-like tail)

• \(c>1\): increasing hazard (wear-out)

Exponentiation parameter \(a\)

Because \(F(x) = F_W(x)^a\):

• \(a>1\) makes the distribution stochastically larger (shifts probability mass to the right).

• \(0<a<1\) makes it stochastically smaller (more mass near 0).

• For integer \(a=m\), \(X\) is the maximum of \(m\) i.i.d. Weibull lifetimes.

Scale and location

• scale rescales the x-axis (units of time/length/etc.).

• loc shifts the support to start at loc.

Below we visualize how the PDF, CDF, and hazard change with parameters.

def plot_pdf_grid(param_sets, *, title: str):
    fig = go.Figure()
    for (a, c, scale) in param_sets:
        rv = exponweib_dist(a, c, scale=scale)
        x_max = rv.ppf(0.995)
        x = np.linspace(1e-6, x_max, 500)
        fig.add_trace(go.Scatter(x=x, y=rv.pdf(x), mode="lines", name=f"a={a}, c={c}, scale={scale}"))

    fig.update_layout(title=title, xaxis_title="x", yaxis_title="pdf")
    fig.show()


def plot_cdf_grid(param_sets, *, title: str):
    fig = go.Figure()
    for (a, c, scale) in param_sets:
        rv = exponweib_dist(a, c, scale=scale)
        x_max = rv.ppf(0.995)
        x = np.linspace(0.0, x_max, 500)
        fig.add_trace(go.Scatter(x=x, y=rv.cdf(x), mode="lines", name=f"a={a}, c={c}, scale={scale}"))

    fig.update_layout(title=title, xaxis_title="x", yaxis_title="cdf")
    fig.show()


def plot_hazard_grid(param_sets, *, title: str):
    fig = go.Figure()
    for (a, c, scale) in param_sets:
        rv = exponweib_dist(a, c, scale=scale)
        x_max = rv.ppf(0.995)
        x = np.linspace(1e-6, x_max, 600)
        h = rv.pdf(x) / rv.sf(x)
        fig.add_trace(go.Scatter(x=x, y=h, mode="lines", name=f"a={a}, c={c}, scale={scale}"))

    fig.update_layout(title=title, xaxis_title="x", yaxis_title="hazard h(x)")
    fig.show()


# Vary 'a' at fixed c
plot_pdf_grid([(0.6, 1.5, 1.0), (1.0, 1.5, 1.0), (2.0, 1.5, 1.0)], title="PDF: effect of a (fixed c=1.5, scale=1)")
plot_cdf_grid([(0.6, 1.5, 1.0), (1.0, 1.5, 1.0), (2.0, 1.5, 1.0)], title="CDF: effect of a (fixed c=1.5, scale=1)")

# Vary 'c' at fixed a
plot_pdf_grid([(1.5, 0.7, 1.0), (1.5, 1.0, 1.0), (1.5, 2.0, 1.0)], title="PDF: effect of c (fixed a=1.5, scale=1)")
plot_hazard_grid([(1.2, 0.7, 1.0), (1.2, 1.0, 1.0), (1.2, 2.0, 1.0)], title="Hazard: effect of c (fixed a=1.2, scale=1)")

# A few parameter combinations that show different hazard shapes
plot_hazard_grid(
    [(0.6, 0.8, 1.0), (1.0, 1.0, 1.0), (2.0, 0.9, 1.0), (2.0, 2.0, 1.0)],
    title="Hazard: different shapes across (a, c)",
)

6) Derivations

We sketch derivations for the expectation/variance and the likelihood.

6.1 Expectation (raw moments)

Work with loc=0, scale=\lambda and consider the \(r\)-th raw moment:

\[\mathbb{E}[X^r] = \int_0^{\infty} x^r\,f(x;a,c,\lambda)\,dx.\]

With

\[f(x) = \frac{a c}{\lambda}\left(\frac{x}{\lambda}\right)^{c-1}\exp\left(-\left(\frac{x}{\lambda}\right)^c\right)\left[1-\exp\left(-\left(\frac{x}{\lambda}\right)^c\right)\right]^{a-1},\]

use the substitution \(y=(x/\lambda)^c\) (so \(x=\lambda y^{1/c}\) and \(dx = (\lambda/c) y^{1/c-1}\,dy\)). The Jacobian cancels the Weibull power term and the integral becomes

\[ \mathbb{E}[X^r] = \lambda^r\,a\int_0^{\infty} y^{r/c}\,e^{-y}\,\bigl(1-e^{-y}\bigr)^{a-1}\,dy. \]

Now expand using the binomial series (valid for \(0<e^{-y}<1\)):

\[ (1-e^{-y})^{a-1} = \sum_{j=0}^{\infty} (-1)^j\binom{a-1}{j}\,e^{-jy}. \]

Swap sum and integral (justified under standard conditions) to get

\[ \mathbb{E}[X^r] = \lambda^r\,a\sum_{j=0}^{\infty} (-1)^j\binom{a-1}{j} \int_0^{\infty} y^{r/c} e^{-(j+1)y}\,dy. \]

Finally, recognize a Gamma integral:

\[ \int_0^{\infty} y^{r/c} e^{-(j+1)y}\,dy = \frac{\Gamma(1+r/c)}{(j+1)^{1+r/c}}, \]

which yields the moment formula used earlier.

6.2 Variance

Compute \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\) from the raw-moment expression and combine:

\[\mathrm{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2.\]

6.3 Likelihood (iid sample)

Given iid data \(x_1,\dots,x_n\) from the standard model with scale \(\lambda\) and shapes \((a,c)\), the likelihood is

\[L(a,c,\lambda) = \prod_{i=1}^n f(x_i; a,c,\lambda).\]

The log-likelihood (often optimized numerically) is

\[ \ell(a,c,\lambda) = n(\log a + \log c - \log \lambda) + (c-1)\sum_{i=1}^n \log(x_i/\lambda) - \sum_{i=1}^n (x_i/\lambda)^c + (a-1)\sum_{i=1}^n \log\bigl(1-\exp(-(x_i/\lambda)^c)\bigr). \]

SciPy’s fit routine maximizes this (with loc/scale included) via numerical optimization.

def exponweib_loglik(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> float:
    """Total log-likelihood using our NumPy logpdf."""
    return float(np.sum(exponweib_logpdf(x, a, c, loc=loc, scale=scale)))


def fit_exponweib_mle_via_minimize(x: np.ndarray, *, loc_fixed: float = 0.0):
    """Simple MLE demo using SciPy optimize on transformed parameters.

    We optimize over (log a, log c, log scale) with loc fixed.
    """
    x = np.asarray(x, dtype=float)

    def nll(theta):
        log_a, log_c, log_scale = theta
        a = np.exp(log_a)
        c = np.exp(log_c)
        scale = np.exp(log_scale)
        return -exponweib_loglik(x, a, c, loc=loc_fixed, scale=scale)

    # crude initial guess from a Weibull fit
    c0, loc0, scale0 = weibull_min.fit(x, floc=loc_fixed)
    theta0 = np.log([1.0, c0, scale0])

    res = optimize.minimize(nll, theta0, method="Nelder-Mead")
    a_hat, c_hat, scale_hat = np.exp(res.x)
    return a_hat, c_hat, loc_fixed, scale_hat, res


# Quick demo on synthetic data
true = dict(a=1.8, c=1.1, loc=0.0, scale=2.5)
x = exponweib_dist.rvs(true["a"], true["c"], loc=true["loc"], scale=true["scale"], size=900, random_state=rng)

(a_hat, c_hat, loc_hat, scale_hat, res) = fit_exponweib_mle_via_minimize(x, loc_fixed=0.0)
print("True params:", true)
print("Minimize MLE:", {"a": a_hat, "c": c_hat, "loc": loc_hat, "scale": scale_hat})
print("Converged:", res.success)

True params: {'a': 1.8, 'c': 1.1, 'loc': 0.0, 'scale': 2.5}
Minimize MLE: {'a': 1.6943268223977832, 'c': 1.1384941757974947, 'loc': 0.0, 'scale': 2.556153488212525}
Converged: True

7) Sampling & Simulation

NumPy-only algorithm (inverse transform)

Because the CDF has a closed-form inverse, sampling is straightforward.

Let \(U\sim\mathrm{Uniform}(0,1)\). Set

\[ X = \text{loc} + \text{scale}\,\left[-\log\left(1-U^{1/a}\right)\right]^{1/c}. \]

Then \(X \sim \mathrm{ExponWeib}(a,c,\text{loc},\text{scale})\).

Implementation detail: for numerical stability,

• compute \(1-\exp(-t)\) as -expm1(-t)

• compute \(\log(1-u)\) as log1p(-u)

• clip \(U\) away from exactly 0 or 1 to avoid log(0).

# Sampling with the NumPy-only sampler

a, c, scale = 1.6, 1.2, 2.0
x = exponweib_rvs_numpy(a, c, scale=scale, size=60_000, rng=rng)

# Compare empirical moments to theory
mean_emp = x.mean()
var_emp = x.var()
mean_theory, var_theory = exponweib_dist.stats(a, c, scale=scale, moments="mv")

print(f"Empirical mean={mean_emp:.4f}, variance={var_emp:.4f}")
print(f"Theory    mean={float(mean_theory):.4f}, variance={float(var_theory):.4f}")

Empirical mean=2.4353, variance=2.7418
Theory    mean=2.4295, variance=2.7334

8) Visualization

We visualize:

• the PDF for different parameter settings

• the CDF (as a sanity check for probability mass)

• Monte Carlo samples vs the theoretical PDF

# PDF / CDF and Monte Carlo comparison

a, c, scale = 1.6, 1.2, 2.0
rv = exponweib_dist(a, c, scale=scale)

x_max = rv.ppf(0.995)
xs = np.linspace(1e-6, x_max, 600)

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=rv.pdf(xs), mode="lines", name="pdf"))
fig.update_layout(title="Exponentiated Weibull PDF", xaxis_title="x", yaxis_title="pdf")
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=rv.cdf(xs), mode="lines", name="cdf"))
fig.update_layout(title="Exponentiated Weibull CDF", xaxis_title="x", yaxis_title="cdf")
fig.show()

# Monte Carlo check
samples = rv.rvs(size=40_000, random_state=rng)
fig = go.Figure()
fig.add_trace(go.Histogram(x=samples, nbinsx=60, histnorm="probability density", name="samples", opacity=0.7))
fig.add_trace(go.Scatter(x=xs, y=rv.pdf(xs), mode="lines", name="theoretical pdf"))
fig.update_layout(title="Samples vs theoretical PDF", xaxis_title="x", yaxis_title="density")
fig.show()

9) SciPy Integration

SciPy exposes this distribution as scipy.stats.exponweib with methods:

• pdf, logpdf, cdf, sf, ppf

• rvs for random variates

• stats / moment / entropy

• fit for maximum-likelihood estimation

A common workflow:

1. Pick a parametric family (here, exponweib).

2. Fit it to data.

3. Check fit visually (histogram, QQ plot) and via likelihood-based criteria.

4. Use the fitted model for inference or simulation.

# SciPy usage: pdf/cdf/rvs/fit

# 1) Create a frozen distribution
rv = exponweib_dist(1.6, 1.2, loc=0.0, scale=2.0)

xs = np.linspace(0, rv.ppf(0.99), 6)
print("x grid:", xs)
print("pdf:", rv.pdf(xs))
print("cdf:", rv.cdf(xs))

# 2) Generate data
x = rv.rvs(size=2_000, random_state=rng)

# 3) Fit parameters (MLE)
# For strictly-positive lifetime data it's often sensible to fix loc=0.
a_hat, c_hat, loc_hat, scale_hat = exponweib_dist.fit(x, floc=0.0)
print("\nFitted params (floc=0):")
print({"a": a_hat, "c": c_hat, "loc": loc_hat, "scale": scale_hat})

# 4) Compare to a Weibull fit (nested model a=1)
a0, c0, loc0, scale0 = exponweib_dist.fit(x, fa=1.0, floc=0.0)
print("\nWeibull-as-exponweib (fa=1) fit:")
print({"a": a0, "c": c0, "loc": loc0, "scale": scale0})

x grid: [0.       1.548109 3.096218 4.644327 6.192436 7.740545]
pdf: [0.       0.295522 0.171121 0.069922 0.024508 0.00785 ]
cdf: [0.       0.351985 0.721415 0.899531 0.967212 0.99    ]

Fitted params (floc=0):
{'a': 1.388222664766818, 'c': 1.2931256583561437, 'loc': 0.0, 'scale': 2.2045665712074545}

Weibull-as-exponweib (fa=1) fit:
{'a': 1.0, 'c': 1.5427975513751417, 'loc': 0.0, 'scale': 2.7025618683212835}

10) Statistical Use Cases

10.1 Hypothesis testing (nested models)

Because the Weibull distribution is the special case \(a=1\), you can test

• \(H_0: a=1\) (Weibull) vs

• \(H_1: a \ne 1\) (exponentiated Weibull)

using a likelihood ratio test (LRT):

\[\Lambda = 2\bigl(\ell(\widehat{a},\widehat{c},\widehat{\lambda}) - \ell(a=1,\widehat{c}_0,\widehat{\lambda}_0)\bigr) \;\approx\; \chi^2_1.\]

10.2 Bayesian modeling

Bayesian modeling treats \((a,c,\lambda)\) as random and combines a prior with the likelihood:

\[p(a,c,\lambda\mid x) \propto p(x\mid a,c,\lambda)\,p(a,c,\lambda).\]

There is no conjugacy here, but generic MCMC (e.g. Metropolis-Hastings) works well.

10.3 Generative modeling

Once fitted (frequentist or Bayesian), the distribution can be used to:

• simulate lifetimes for stress testing

• generate synthetic positive-valued data with realistic skew/tails

• build parametric simulators inside larger pipelines

# 10.1 Likelihood ratio test: Weibull (a=1) vs Exponentiated Weibull

# Simulated data under the alternative
true = dict(a=1.8, c=1.1, loc=0.0, scale=2.5)
x = exponweib_dist.rvs(true["a"], true["c"], loc=true["loc"], scale=true["scale"], size=900, random_state=rng)

# Fit alternative (a free) with loc fixed to 0
(a1, c1, loc1, scale1) = exponweib_dist.fit(x, floc=0.0)
ll1 = float(np.sum(exponweib_dist.logpdf(x, a1, c1, loc=0.0, scale=scale1)))

# Fit null (a=1) with loc fixed to 0
(a0, c0, loc0, scale0) = exponweib_dist.fit(x, fa=1.0, floc=0.0)
ll0 = float(np.sum(exponweib_dist.logpdf(x, a0, c0, loc=0.0, scale=scale0)))

lr_stat = 2.0 * (ll1 - ll0)
p_value = chi2.sf(lr_stat, df=1)

print("Alt fit:", {"a": a1, "c": c1, "scale": scale1})
print("Null fit:", {"a": a0, "c": c0, "scale": scale0})
print(f"LR statistic={lr_stat:.3f}, p-value={p_value:.4g}")

Alt fit: {'a': 2.2710907010507073, 'c': 0.9742363049697826, 'scale': 2.082356652940917}
Null fit: {'a': 1.0, 'c': 1.4663239585093428, 'scale': 3.7688940378192104}
LR statistic=21.856, p-value=2.94e-06

# 10.2 Bayesian modeling (toy example): random-walk Metropolis on log-parameters

x = exponweib_dist.rvs(1.8, 1.1, scale=2.5, size=500, random_state=rng)

# Prior: independent normals on log-parameters (broad, weakly informative)
prior_mu = np.array([0.0, 0.0, 0.0])
prior_sigma = np.array([1.5, 1.5, 1.5])


def log_prior(theta: np.ndarray) -> float:
    d = theta.size
    return float(
        -0.5 * np.sum(((theta - prior_mu) / prior_sigma) ** 2)
        - np.sum(np.log(prior_sigma))
        - 0.5 * d * np.log(2 * np.pi)
    )


def log_lik(theta: np.ndarray) -> float:
    log_a, log_c, log_scale = theta
    a, c, scale = np.exp([log_a, log_c, log_scale])
    return float(np.sum(exponweib_dist.logpdf(x, a, c, loc=0.0, scale=scale)))


def log_post(theta: np.ndarray) -> float:
    return log_lik(theta) + log_prior(theta)


n_steps = 6_000
burn = 1_500
step_scale = np.array([0.08, 0.08, 0.08])

# Initialize at the MLE (good starting point)
a_mle, c_mle, loc_mle, scale_mle = exponweib_dist.fit(x, floc=0.0)
cur = np.log([a_mle, c_mle, scale_mle])
cur_lp = log_post(cur)

chain = np.empty((n_steps, 3), dtype=float)
accepted = 0

for t in range(n_steps):
    prop = cur + step_scale * rng.normal(size=3)
    prop_lp = log_post(prop)

    if np.log(rng.uniform()) < (prop_lp - cur_lp):
        cur, cur_lp = prop, prop_lp
        accepted += 1

    chain[t] = cur

acc_rate = accepted / n_steps
print(f"Acceptance rate: {acc_rate:.3f}")

post = chain[burn:]
a_s, c_s, scale_s = np.exp(post.T)

summary = {
    "a_mean": float(a_s.mean()),
    "c_mean": float(c_s.mean()),
    "scale_mean": float(scale_s.mean()),
    "a_95%": tuple(np.quantile(a_s, [0.025, 0.975])),
    "c_95%": tuple(np.quantile(c_s, [0.025, 0.975])),
    "scale_95%": tuple(np.quantile(scale_s, [0.025, 0.975])),
}
summary

Acceptance rate: 0.178

{'a_mean': 1.4587089908602129,
 'c_mean': 1.2533464400252319,
 'scale_mean': 3.084652557279298,
 'a_95%': (1.0627486616378665, 2.2450850386814643),
 'c_95%': (1.0083703414202325, 1.4947448812498432),
 'scale_95%': (2.1649431209050984, 3.8530214219000483)}

# Simple posterior trace plots

fig = go.Figure()
fig.add_trace(go.Scatter(y=a_s, mode="lines", name="a"))
fig.update_layout(title="Posterior trace: a", xaxis_title="iteration", yaxis_title="a")
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(y=c_s, mode="lines", name="c"))
fig.update_layout(title="Posterior trace: c", xaxis_title="iteration", yaxis_title="c")
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(y=scale_s, mode="lines", name="scale"))
fig.update_layout(title="Posterior trace: scale", xaxis_title="iteration", yaxis_title="scale")
fig.show()

# 10.3 Generative modeling: posterior predictive vs fitted MLE

# Use the same x from the Bayesian section

# MLE predictive
rv_mle = exponweib_dist(a_mle, c_mle, loc=0.0, scale=scale_mle)
x_mle = rv_mle.rvs(size=15_000, random_state=rng)

# Posterior predictive (draw parameters, then draw a sample)
idx = rng.integers(0, len(a_s), size=400)
pp_samples = []
for i in idx:
    rv_i = exponweib_dist(a_s[i], c_s[i], loc=0.0, scale=scale_s[i])
    pp_samples.append(rv_i.rvs(size=40, random_state=rng))
pp_samples = np.concatenate(pp_samples)

fig = go.Figure()
fig.add_trace(go.Histogram(x=x, nbinsx=60, histnorm="probability density", name="observed", opacity=0.6))
fig.add_trace(go.Histogram(x=x_mle, nbinsx=60, histnorm="probability density", name="MLE predictive", opacity=0.5))
fig.add_trace(go.Histogram(x=pp_samples, nbinsx=60, histnorm="probability density", name="Posterior predictive", opacity=0.5))
fig.update_layout(title="Observed vs predictive distributions", barmode="overlay", xaxis_title="x", yaxis_title="density")
fig.show()

11) Pitfalls

• Invalid parameters: require \(a>0\), \(c>0\), scale>0.

• Behavior at \(x\approx 0\): depending on \(c\) and \(a\), the density may go to 0 or blow up; use logpdf for stability.

• CDF cancellation: for very small \(x\), \(1-\exp(-x^c)\) suffers cancellation; use expm1 (-expm1(-t)).

• Inverse-CDF sampling: avoid exact 0/1 uniforms; use clip and log1p.

• Fitting: unconstrained loc can drift negative even for strictly-positive data; fix floc=0 when appropriate.

• Moment calculations: the binomial series can require many terms for extreme parameters; prefer SciPy numerical moments or Monte Carlo when in doubt.

• MGF existence: for \(c\le 1\), the MGF may not exist for positive \(t\); use the characteristic function or Laplace transform instead.

12) Summary

• exponweib is a continuous distribution on \([0,\infty)\) (standard form) with shape parameters \(a>0\) and \(c>0\).

• It is defined by exponentiating a Weibull CDF: \(F(x) = [1-\exp(-x^c)]^a\) (with optional loc/scale).

• It generalizes Weibull (a=1) and can represent diverse hazard shapes important in reliability/survival analysis.

• Raw moments admit a useful Gamma + binomial series representation; SciPy provides numerical stats and entropy.

• Sampling is easy via inverse transform and can be implemented with NumPy only.

References

• SciPy documentation: scipy.stats.exponweib.

• Wikipedia: “Exponentiated Weibull distribution”.